# =============================================================== #
# Reproducibility of:
# Drivers and Trends for the Equality of Opportunity for Sexual
# and Gender Minorities: A Panel Approach (EQOSOGI 2024)
#
# Script: Cluster features tables (one table per cluster)
# =============================================================== #

library(readxl)
library(tidyverse)
library(cluster)
library(e1071)
library(clusterCrit)
library(writexl)

# Rutas
project_root    <- getwd()
intermediate_dir <- file.path(project_root, "intermediate files")
outputs_dir      <- file.path(project_root, "outputs")

# ------------------------------------------------------ #
# 1. Cargar datos en panel y construir features por país #
# ------------------------------------------------------ #

data_long <- read_excel(file.path(intermediate_dir, "eqosogi_score.xlsx"))

features_df <- data_long %>%
  group_by(country) %>%
  arrange(year) %>%
  nest() %>%
  mutate(
    slope_total = map_dbl(data, ~ coef(lm(es_ ~ year, data = .x))[2]),
    delta       = map_dbl(data, ~ last(.x$es_) - first(.x$es_)),
    last_value  = map_dbl(data, ~ last(.x$es_)),
    sd          = map_dbl(data, ~ sd(.x$es_)),
    mean        = map_dbl(data, ~ mean(.x$es_)),
    max_jump    = map_dbl(data, ~ max(diff(.x$es_))),
    slope_late  = map_dbl(data, ~ coef(lm(es_ ~ year,
                                          data = dplyr::filter(.x, year >= 2001)))[2]),
    delta_late  = map_dbl(data, ~ .x$es_[.x$year == 2024] -
                            .x$es_[.x$year == 2001])
  ) %>%
  select(-data) %>%
  ungroup()

# -------------------------------------- #
# 2. k-means y asignación de 3 clusters  #
# -------------------------------------- #

features_scaled <- scale(select(features_df, -country))

set.seed(123)
kmeans_result <- kmeans(features_scaled, centers = 3)

# Guardamos cluster numérico y lo re-etiquetamos como en el script original
features_df <- features_df %>%
  mutate(
    cluster_num = kmeans_result$cluster,
    cluster = case_when(
      cluster_num == 1 ~ "Cluster 1",  # flat / limited progress
      cluster_num == 2 ~ "Cluster 3",  # steepest growth
      cluster_num == 3 ~ "Cluster 2"   # intermediate
    ),
    cluster = factor(cluster, levels = c("Cluster 1", "Cluster 2", "Cluster 3"))
  )

# ----------------------------------------------- #
# 3. Construir tablas de features por cada cluster #
# ----------------------------------------------- #

feature_cols <- c("slope_total", "delta", "last_value",
                  "sd", "mean", "max_jump", "slope_late", "delta_late")

make_cluster_table <- function(df, cluster_label) {
  # Filtrar cluster y ordenar por país
  df_cluster <- df %>%
    filter(cluster == cluster_label) %>%
    arrange(country) %>%
    select(Country = country, all_of(feature_cols))
  
  # Fila de promedio
  avg_row <- df_cluster %>%
    summarise(across(all_of(feature_cols), ~ mean(.x, na.rm = TRUE))) %>%
    mutate(Country = "Average") %>%
    select(Country, all_of(feature_cols))
  
  # Unir países + promedio y redondear a 4 decimales
  bind_rows(df_cluster, avg_row) %>%
    mutate(across(all_of(feature_cols), ~ round(.x, 4)))
}

table_cluster1 <- make_cluster_table(features_df, "Cluster 1")
table_cluster2 <- make_cluster_table(features_df, "Cluster 2")
table_cluster3 <- make_cluster_table(features_df, "Cluster 3")

# -------------------------------------------- #
# 4. Exportar a Excel (un sheet por cluster)   #
# -------------------------------------------- #

output_path <- file.path(outputs_dir, "Tables_A2-1_to_A2-3.xlsx")

write_xlsx(
  list(
    "Cluster 1" = table_cluster1,
    "Cluster 2" = table_cluster2,
    "Cluster 3" = table_cluster3
  ),
  path = output_path
)
